library(tidyverse)
library(ggplot2)
library(lme4)
library(dotwhisker)
library(broom)
library(ggrepel)
library(stargazer)

# carregando dados
cleadefr <- read.csv("clea_de-fr.csv")
essdefr <- read.csv("ess_de_fr_12.csv")
de17 <- read.csv("2017.csv")
bwl17 <- read.csv("btw17_wbz_zweitstimmen.csv", sep = ";")
bwl13 <- read.csv("BTW13_Zweitstimmen_Wahlbezirke.csv", sep = ";")


# renomeando observações
cleadefr$pty_n[cleadefr$pty_n == "nationaldemokratische partei deutschlands"] <- "NPD"
cleadefr$pty_n[cleadefr$pty_n == "National Front"] <- "front national"

##### descriptive analysis #####
attach(cleadefr)

# line plot 1 - FN average vote share by election year 
line1 <- ggplot(cleadefr[cleadefr$pty_n == "front national",], aes(x = yr, y = pvs1)) +
  theme_classic() +
  scale_x_continuous(limits = c(1978, 2017), breaks = c(1978, 1981, 1986, 1988,
                                                        1993, 1997, 2002, 2007, 
                                                        2012, 2017)) +
  theme(plot.title = element_text(face = "bold", hjust = .5, size = 12),
        axis.text.x = element_text(angle = 90)) +
  labs(x = "Ano", y = "Média de votos (%)") +
  ggtitle("Média de votos do FN por ano de eleição") +
  stat_summary(fun.y = mean, geom = "line")

ggsave("FN_vote_timeseries.png", plot = line1, device = "png", width = 8, height = 5)

# line plot 2 - constituency level vote share by year of election
line2 <- ggplot(cleadefr[cleadefr$pty_n == "front national",], aes(x = yr, y = pvs1)) +
  geom_line(aes(group = cst), alpha = .1) +
  geom_smooth(method = "lm", colour = "red", se = F) +
  theme_classic() +
  ggtitle("Porcentagem de votos do FN nas constituências por ano de eleição") +
  labs(x = "Ano", y = "Porcentagem de votos") +
  theme(plot.title = element_text(size = 12, hjust = .5, face = "bold")) +
  scale_x_continuous(limits = c(1978, 2017), breaks = c(1978, 1981, 1986, 1988,
                                                        1993, 1997, 2002, 2007, 
                                                        2012, 2017))

ggsave("FN_cst_level_timeseries.png", plot = line2, device = "png", width = 8, height = 5)

# line plot 3 - NPD average vote share by year
line3 <- ggplot(cleadefr[cleadefr$ctr_n == "Germany" & cleadefr$pty_n %in% c("nationaldemokratische partei deutschlands",
                                                           "NPD"),], 
                aes(x = yr, y = pvs1)) +
  theme_classic() +
  stat_summary(fun.y = mean, geom = "line") +
  scale_x_continuous(limits = c(1965, 2013), breaks = c(1965, 1969, 1972, 1976, 1980,
                                                        1983, 1987, 1990, 1994, 1998,
                                                        2002, 2005, 2009, 2013)) +
  theme(plot.title = element_text(face = "bold", size = 12, hjust = .5)) +
  ggtitle("Média de votos do NPD por ano de eleição") +
  labs(x = "Ano", y = "Média de votos (%)")

ggsave("NPD_timeseries.png", plot = line3, device = "png", width = 8, height = 5)

# line plot 4 - NPD constituency level vote share by year of election
line4 <- ggplot(cleadefr[cleadefr$ctr_n == "Germany" & cleadefr$pty_n %in% c("nationaldemokratische partei deutschlands",
                                                           "NPD"),], 
                aes(x = yr, y = pvs1)) +
  geom_line(aes(group = cst), alpha = .1) +
  geom_smooth(method = "lm", colour = "red", se = F, na.rm = T) +
  theme_classic() +
  ggtitle("Porcentagem de votos do NPD nas constituências por ano de eleição") +
  labs(x = "Ano", y = "Porcentagem de votos") +
  theme(plot.title = element_text(size = 12, hjust = .5, face = "bold")) +
  scale_x_continuous(limits = c(1965, 2013), breaks = c(1965, 1969, 1972, 1976, 1980,
                                                        1983, 1987, 1990, 1994, 1998,
                                                        2002, 2005, 2009, 2013))

ggsave("NPD_cst_level_timeseries.png", plot = line4, device = "png", width = 8, height = 5)

# boxplot 1 - AfD vote distribution 2013
cleadefr$east <- as.factor(ifelse(cleadefr$sub %in% c("Thüringen", "Sachsen-Anhalt",
                                                      "Brandenburg", "Sachsen", 
                                                      "Mecklenburg-Vorpommern"), 1, 0))

box1 <- ggplot(cleadefr[ctr_n == "Germany" & pty_n %in% c("AfD", "alternative for germany"),], 
               aes(y = pvs1, x = factor(sub), fill = east)) +
  geom_boxplot(width = .5) +
  coord_flip() +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 12, hjust = .5),
        legend.title = element_text(face = "bold", size = 10, hjust = .5)) +
  ggtitle("Distribuição dos votos na AfD por região nas eleições de 2013") +
  labs(x = "Região", y = "Média de votos (%)") +
  scale_color_continuous(name = "Leste")

ggsave("boxplot_AfD_2013.png", plot = box1, device = "png", width = 8, height = 5)

# boxplot 2 - AfD vote distribution 2017
box2 <- ggplot(de17[de17$Registered.Electors == "Alternative for Germany (AfD)",],
               aes(y = X, x = factor(X2017))) +
  geom_boxplot(width = .3) +
  geom_point(position = position_jitter(.15)) +
  stat_summary(fun.y = mean, geom = "point", colour = "red", size = 1.5) +
  theme_minimal() +
  ggtitle("Boxplot da distribuição dos votos da AfD nas eleições de 2017") +
  theme(plot.title = element_text(size = 12, hjust = .5, face = "bold")) +
  labs(x = "Eleição", y = "Média de votos (%)") +
  geom_label_repel(data = subset(de17[de17$Registered.Electors == "Alternative for Germany (AfD)",], 
                                 X > 22), 
                   aes(factor(X2017), X, label = DE),
                   box.padding   = 0.35, 
                   point.padding = 0.5,
                   segment.color = 'grey50', direction = "x")

ggsave("boxplot_AfD_2017_general.png", plot = box2, device = "png", width = 8, height = 5)  

# boxplot 3 - AfD vote distribution 2017 by land

# creating variable for east/west land
bwl17$east <- as.factor(ifelse(bwl17$Land %in% c(13, 12, 15, 14, 16), 1, 0))

# boxplot - booyah, motherfuckeeeer
box3 <- ggplot(bwl17, aes(x = factor(Land), y = AfD, fill = east)) +
  geom_boxplot() +
  theme_minimal() +
  theme(plot.title = element_text(size = 12, face = "bold", hjust = .5)) +
  ggtitle("Distribuição dos votos na AfD por Land na Alemanha") +
  labs(x = "Região", y = "Número de votos AfD")

ggsave("boxplot_AfD_2017.png", plot = box3, device = "png", width = 8, height = 5)

# função para nomear outliers
is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

bwl17 %>% mutate(outlier = ifelse(is_outlier(AfD), AfD, as.numeric(NA))) %>%
ggplot(., aes(x = factor(Land), y = AfD)) +
  geom_boxplot() +
  geom_text(aes(label = outlier))

detach(cleadefr)

##### modelos #####
attach(essdefr)

# variáveis demográficas
logit1 <- glmer(vote_rrp ~ gndr + agea + (1 | cntry), family = "binomial")
summary(logit1)

# variáveis socioeconômicas
logit2 <- glmer(vote_rrp ~ eduyrs + blue_c + uemp5yr + (1 | cntry), family = "binomial")
summary(logit2)

# modelo geral com controles
logit3 <- glmer(vote_rrp ~ gndr + agea + eduyrs + blue_c + uemp5yr + lrscale +
                  (1 | cntry), family = "binomial")
summary(logit3, corr = F)

# modelos por país: Alemanha
logit4 <- glmer(data = essdefr[essdefr$cntry == "DE",], vote_rrp ~ gndr + agea +
                  eduyrs + blue_c + uemp5yr + lrscale + (1 | intewde), 
                family = "binomial")
summary(logit4, corr = F)

# Alemanha não multinível
logit7 <- glm(data = essdefr[essdefr$cntry == "DE",], vote_rrp ~ gndr + agea +
                eduyrs + blue_c + uemp5yr + lrscale, family = "binomial")
summary(logit7)

# França
logit5 <- glm(data = essdefr[essdefr$cntry == "FR",], vote_rrp ~ gndr + agea +
                  eduyrs + blue_c + uemp5yr + lrscale, 
                family = "binomial")
summary(logit5, corr = F)

# França não multinível
logit6 <- glm(data = essdefr[essdefr$cntry == "FR",], vote_rrp ~ gndr + agea +
                eduyrs + blue_c + uemp5yr + lrscale, family = "binomial")
summary(logit6)

##### plot dos coeficientes #####
# renomeando modelos multinível
modelos1 <- list(logit1, logit2, logit3, logit4)

names(modelos1) <- c("Demográficas", "Socioeconômicas", 
                    "Geral", "Alemanha")

# renomeando modelos logísticos
modelos2 <- list(logit6, logit7)
names(modelos2) <- c("França", "Alemanha")

# tentativa de renomear usando o broom e o dplyr 
m1_tidy <- tidy(logit1)
m2_tidy <- tidy(logit2)
m3_tidy <- tidy(logit3)
m4_tidy <- tidy(logit4)

# modificando nomes dos modelos multinível
m1_tidy <-  m1_tidy %>% mutate(model = "(1) Demográficas")
m2_tidy <- m2_tidy %>% mutate(model = "(2) Socieconômicas")
m3_tidy <- m3_tidy %>% mutate(model = "(3) Geral")
m4_tidy <- m4_tidy %>% mutate(model = "(4) Alemanha")

# binding
allmodels1 <- bind_rows(m1_tidy, m2_tidy, m3_tidy, m4_tidy)

# modificando nomes dos modelos
m6_tidy <- tidy(logit6)
m7_tidy <- tidy(logit7)

# modificando nomes
m6_tidy <- m6_tidy %>% mutate(model = "França")
m7_tidy <- m7_tidy %>% mutate(model = "Alemanha")

# binding
allmodels2 <- bind_rows(m6_tidy, m7_tidy)

# renomeando variáveis
allmodels1 <- allmodels1 %>%
  relabel_predictors(c("sd_(Intercept).cntry" = "s.d. País",
                       "sd_(Intercept).intewde" = "s.d. Leste/Oeste Alemanha",
                       "sd_(Intercept).regionfr" = "s.d. Região França",
                       "gndr" = "Gênero",
                       "agea" = "Idade (num.)",
                       "eduyrs" = "Educação (anos)",
                       "blue_c" = "Blue-collar",
                       "uemp5yr" = "Desempregado (5 anos)",
                       "lrscale" = "Esquerda-direita"))

# renomeando variáveis 
allmodels2 <- allmodels2 %>% 
  relabel_predictors(c("gndr" = "Gênero",
                       "agea" = "Idade (anos)",
                       "eduyrs" = "Educação (anos)",
                       "blue_c" = "Blue-collar",
                       "uemp5yr" = "Desempregado (5 anos)",
                       "lrscale" = "Esquerda-direita"))

# modelos multinível
coef1 <- dwplot(modelos1, show_intercept = F, 
              vline = geom_vline(xintercept = 0, colour = "grey 60", 
                                 linetype = 2)) +
  ggtitle("Coeficientes dos Modelos multinível") +
  labs(x = "Coeficientes com intervalo de confiança a 95%") +
  theme_classic() +
  theme(plot.title = element_text(size = 12, face = "bold", hjust = .5),
        legend.title = element_text(size = 10, face = "bold"),
        legend.position = "right", legend.text.align = .5) +
  scale_color_discrete("Variáveis")
  
# modelos não multinível
coef2 <- dwplot(modelos2, show_intercept = F, 
                vline = geom_vline(xintercept = 0, colour = "grey 60", 
                                   linetype = 2)) +
  ggtitle("Coeficientes dos Modelos logísticos") +
  labs(x = "Coeficientes com intervalo de confiança a 95%") +
  theme_classic() +
  theme(plot.title = element_text(size = 12, face = "bold", hjust = .5),
        legend.title = element_text(size = 10, face = "bold"),
        legend.position = "right", legend.text.align = .5) +
  scale_color_discrete("Países")

# modelos multinível tidy
coef3 <- dwplot(allmodels1, show_intercept = F, 
                vline = geom_vline(xintercept = 0, colour = "grey 60", 
                                   linetype = 2)) +
  ggtitle("Coeficientes dos Modelos multinível") +
  labs(x = "Coeficientes com intervalo de confiança a 95%") +
  theme_classic() +
  theme(plot.title = element_text(size = 12, face = "bold", hjust = .5),
        legend.title = element_text(size = 10, face = "bold"),
        legend.position = "right", legend.text.align = .5) +
  scale_color_discrete("Variáveis")

ggsave("coef_multilevel.png", plot = coef3, device = "png", width = 6, height = 5)

# modelos normais tidy
coef4 <- dwplot(allmodels2, show_intercept = F, 
                vline = geom_vline(xintercept = 0, colour = "grey 60", 
                                   linetype = 2)) +
  ggtitle("Coeficientes dos Modelos logísticos") +
  labs(x = "Coeficientes com intervalo de confiança a 95%") +
  theme_classic() +
  theme(plot.title = element_text(size = 12, face = "bold", hjust = .5),
        legend.title = element_text(size = 10, face = "bold"),
        legend.position = "right", legend.text.align = .5) +
  scale_color_discrete("Países")

ggsave("coef_gen.png", plot = coef4, device = "png", width = 6, height = 5)

##### tabela com os coeficientes #####
stargazer(logit1, logit2, logit3, logit4,
          title = "Resultados - modelos logit multinível", 
          dep.var.labels = "Voto em partidos populistas",
          column.labels = c("Demográficas", "Socioeconômicas", "Geral", "Alemanha"),
          covariate.labels = c("Gênero", "Idade", "Educação", "Blue-collar", 
                               "Desempregado", "Esquerda-direita", "Constante"),
          model.numbers = T,
          mean.sd = T, min.max = T, median = T, iqr = T,
          type = "html", align = T, out = "results1.htm")

stargazer(logit6, logit7,
          title = "Resultados - modelos logit", 
          dep.var.labels = "Voto em partidos populistas",
          column.labels = c("França", "Alemanha"),
          covariate.labels = c("Gênero", "Idade", "Educação", "Blue-collar", 
                               "Desempregado", "Esquerda-direita", "Constante"),
          model.numbers = T,
          mean.sd = T, min.max = T, median = T, iqr = T,
          type = "html", align = T, out = "results2.htm")


##### salvando plots #####
setwd("C:/Users/test/Desktop/Seminário CPRI 2019/plots")

ggsave()